{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn import datasets\n",
    "boston = datasets.load_boston()\n",
    "X = boston['data']\n",
    "y = boston['target']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `BayesianRegression` class estimates the regression coefficients using \n",
    "\n",
    "$$\n",
    "\\left(\\frac{1}{\\sigma^2}\\bX^\\top\\bX + \\frac{1}{\\tau} I\\right)^{-1}\\frac{1}{\\sigma^2}\\bX^\\top\\by.\n",
    "$$\n",
    "\n",
    "Note that this assumes $\\sigma^2$ and $\\tau$ are known. We can determine the influence of the prior distribution by manipulationg $\\tau$, though there are principled ways to choose $\\tau$. There are also principled Bayesian methods to model $\\sigma^2$ (see [here](https://www.statlect.com/fundamentals-of-statistics/Bayesian-regression)), though for simplicity we will estimate it with the typical OLS estimate:\n",
    "\n",
    "$$\n",
    "\\hat{\\sigma}^2 = \\frac{SSE}{N - (D + 1)},\n",
    "$$\n",
    "\n",
    "where $SSE$ is the sum of squared errors from an ordinary linear regression, $N$ is the number of observations, and $D$ is the number of predictors. Using the linear regression model from chapter 1, this comes out to about 11.8."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BayesianRegression:\n",
    "    \n",
    "    def fit(self, X, y, sigma_squared, tau, add_intercept = True):\n",
    "        \n",
    "        # record info\n",
    "        if add_intercept:\n",
    "            ones = np.ones(len(X)).reshape((len(X),1))\n",
    "            X = np.append(ones, np.array(X), axis = 1)\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "        \n",
    "        # fit\n",
    "        XtX = np.dot(X.T, X)/sigma_squared\n",
    "        I = np.eye(X.shape[1])/tau\n",
    "        inverse = np.linalg.inv(XtX + I)\n",
    "        Xty = np.dot(X.T, y)/sigma_squared\n",
    "        self.beta_hats = np.dot(inverse , Xty)\n",
    "        \n",
    "        # fitted values\n",
    "        self.y_hat = np.dot(X, self.beta_hats)\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's fit a Bayesian regression model on the {doc}`Boston housing </content/appendix/data>` dataset. We'll use $\\sigma^2 = 11.8$ and $\\tau = 10$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_squared = 11.8\n",
    "tau = 10\n",
    "model = BayesianRegression()\n",
    "model.fit(X, y, sigma_squared, tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below plot shows the estimated coefficients for varying levels of $\\tau$. A lower value of $\\tau$ indicates a stronger prior, and therefore a greater pull of the coefficients towards their expected value (in this case, 0). As expected, the estimates approach 0 as $\\tau$ decreases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABJQAAAE0CAYAAAB6jmkMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZgtdXkn8O8LiKhg1OAGcl3irmMQb4xGA2rcMHFN3GJMJjESjQwxjkscnZjNuIwxjjNGhWjMOGp0VNxQUB8FlzEuKCC4ZBBFtkRRUUCEiO/8cc7Vpu0+fere7ntOn/58nqee21V1uur9naVe+HZVneruAAAAAMC09ph1AQAAAABsLgIlAAAAAAYRKAEAAAAwiEAJAAAAgEEESgAAAAAMIlACAAAAYBCBEgAAAACDCJQAAAAAGESgxIapqjOq6l6zrmOIqrpNVX2uqi6uqqNWmJ9qTJtp7GvVWlVfq6r77saSgC1qMx07d9A3VlyvbwC7xWY6du6gb6y4Xt/YpARKc2z8wbqsqi6pqn+tqtdV1b6zrmta3X2H7j5xI7ZdVb9ZVZ8ZPzcXVNX7quqe67DpZyY5sbv36+6XL5+fdkzrOfaNPsAurXWzHMyr6sjx6395Vb1uhfXXq6pjq+rSqjq7qn5zyHrYrPSN1ekb62cr9g1YVPrG6vSN9bOIfYMRgdL8e3B375vk4CR3TvLs9d5BVe213tvcSFX1tCQvS/LXSW6YZFuSv0vy0HXY/E2TnDFhnvlwfpK/SvLaVda/IskVGb0/HpfklVV1hwHrYTPTN5bRN8iu9w1YZPrGMvoGWbtvkCTdbZrTKcnXktx3yfyLkxy3ZP6AJG9L8s0kX01y1JJ1hyT5XJKLk/yfJG9O8lfLtv2sJKcluTzJXmts71lJzhtv78tJfmXS8uX1J7ldkhOTXJTRAfMhyx739HEt3x3Xus8qz8nPJLkkySMnPG+T9jVpjB9KcmWSH4z3sXz+1iu8Jgcleft4e99K8j9Xee0m7XfV8Sd5fZIfJblsXMMzJz3nS7b5u0nevWT+zCRvWTJ/TpKDl9a60r6Gvj7jx39x/PtXjKdLxtPtNuAz8ldJXrds2bXG+731kmWvT/LCadabTJt5WuHYo2/oG/rGVfc1uG+YTIs8rXDs0Tf0DX3jqvv6qb5hWvL8zLoA04QX56oHyJsk+XyS/z6e3yPJyUn+NMneSW6R5KwkDxjPn53kj5JcLckjxh+05Qf4U8YHqGussb3bjA8IB4x/92ZJfm615cvrH9dwZpL/Mt72fcYHp9ssedynxgfB640PEE9a5Tl5YJIfJtlrlfWr7mvSGJf8/olJfn/C/NLXZM8kpyb524z+Y3SfJPdc4XET97vW+Jdta+JzvuR3bpFRg9sjyY3H74fzlqz7TpI9Vtj+j39etv+pXp9lv/eaJM9ZZd17xvWtNL1nwGdkpf8xuHOSy5Yte3rGDW+t9SbTZp6ib6z0nOgb+sbS7QzuGybTIk/RN1Z6TvQNfWPpdgRKEyaXvM2/d1TVxRl9qL+R5Hnj5b+Q5Prd/RfdfUV3n5XkmCSPSXK3jP4C8PLu/vfufntGH9DlXt7d53T3ZWts78okV09y+6q6Wnd/rbu/MmH5cndLsm9Gf+m7ors/lNEH/LHLajm/u7+d5N0ZnXK7kp9NcmF3/3CV9ZP2NWmMO+OuGR30ntHdl3b3D7r7Yys8bpr9Tjv+qZ7z8T4uHm/nsCQnJDmvqm47nv9od/9owFinrW+pOyU5faUV3f1r3X2dVaZfG1DXSvbN6C8bS303yX5TrofNTt+4Kn1D31iLvsBWp29clb6hbzClTXUt6xb1sO7+YFUdluSNSfbPKFW9aZIDquqiJY/dM8lHMzronNc9ilTHzllh20uXrbq97j6zqp6a5M+S3KGqTkjytAnLz1+2nwOSnLPsgHJ2kgOXzP/rkp+/P/6dlXwryf5VtdcqB/lJ+5r0nO2Mg5KcPaHZ7DDNfqca/4DnPElOSnKvJLcc/3xRRgf3u4/nh5j29UmSVNUeSW6fVQ7wG+ySJNdetuzaGTW8adbDZqdvXJW+oW+sRV9gq9M3rkrf0DeYkjOUNonuPinJ65K8ZLzonCRfXZa07tfdD0pyQZIDq6qWbOKglTa75OdJ20t3v7G775nRwaqTvGjS8mXOT3LQ+EO/w7aMrssd6hMZXWP8sFXWT9rXxDHuhHOSbJviJoO7ut++ysx0z3nykwP8L49/PimjA/xhWf0A36ssH2pbRseXs1ZaOf6WjEtWmd63i/v+lyR7VdWtliz7+fzkZodrrYeFoG/8mL6hb6xFX4DoG0voG/oGUxIobS4vS3K/qjo4o1NKv1dVz6qqa1TVnlV1x6r6hYwOglcmObKq9qqqh2Z0uuQkq26vqm5TVfepqqtndHC9LMmVqy1fYdufTHJpkmdW1dWq6l5JHpzkn4Y+Ad393YyuDX5FVT2sqq453ubhVfXiNfY16TnbGZ/KqJm+sKquVVX7VNU9Vnncruz33zK6DjkDnvNkdBC/d5JrdPe5Gf2F4oEZncb7ubX2tYuundHrsPdKK7v78O7ed5Xp8LU2Pn5f75PRX172HD/3e423fWlGNy78i/Hrco+MvpHj9dOshwWjb+gb+kZ2rW/AFqNv6Bv6Rib3DX5CoLSJdPc3k/yvJP+1u6/M6MB1cEZ38b8wyd8n+ZnuviKjG+M9IaPTDn8ro+t6L5+w7VW3l9E1tC8cL/vXJDfI6CZ0qy1fvu0rkjwkyeHjx/5dkt/u7i/t5PPw0iRPS/LcjL7F4JwkRyZ5x6R9rTHGnaljx/ZumeTrSc5N8ugJj9vZ/b4gyXNrdArrozPFcz7e779kdBr/R8fz38sowf/4uKaJ+6qqp09Z30q+mNENBL9To+uo19tzM2puf5LR+/uy8bId/jCjmz9+I8mbkjy5u88YsB4Wgr7x4+3pG/rGrvYN2BL0jR9vT9/QN9bqGySp7vU644x5VlWfTPKq7v6HWdcCwPzTNwAYQt+ArccZSguqqg6rqhuNT9X7nYzugH/8rOsCYD7pGwAMoW8ArgFcXLdJ8paMvtLyK0l+o7svmG1JAMwxfQOAIfQN2OJc8gYAAADAIC55AwAAAGAQgRIAAAAAgyzEPZQe+MAH9vHHu/8bwCpq1gXMmj4BsKot3yMSfQJgglX7xEKcoXThhRfOugQA5pg+AcAk+gTAcAsRKAEAAACw+wiUAAAAABhEoAQAAADAIAIlAAAAAAYRKAEAAAAwiEAJAAAAgEEESgAAAAAMIlACAAAAYBCBEgAAAACD7DXrAgAAAGCWjjnus+u6vSf+6iHruj2YR3MbKFXV15JcnOTKJD/s7u2zrQgAAACAZI4DpbF7d/eFsy4CAAAAgJ9wDyUAAAAABpnnQKmTvL+qTq6qI2ZdDAAAAAAj83zJ2z26+/yqukGSD1TVl7r7IztWjkOmI5Jk27Zts6oRgDmlTwAwiT4BsGvmNlDq7vPH/36jqo5NctckH1my/ugkRyfJ9u3beyZFAjC39AmW8u09wHL6BMCumctL3qrqWlW1346fk9w/yemzrQoAAACAZH7PULphkmOrKhnV+MbuPn62JQEAAACQzGmg1N1nJfn5WdcBAAAAwE+by0veAAAAAJhfAiUAAAAABhEoAQAAADCIQAkAAACAQQRKAAAAAAwyl9/yBkMdc9xn13V7T/zVQ9Z1ewAAALBInKEEAAAAwCACJQAAAAAGESgBAAAAMIhACQAAAIBBBEoAAAAADCJQAgAAAGAQgRIAAAAAgwiUAAAAABhEoAQAAADAIAIlAAAAAAYRKAEAAAAwiEAJAAAAgEEESgAAAAAMIlACAAAAYJC9Zl0As/XJz35tXbf3i4fcbF23BwAAAMwfgRIAsKL1/qND4g8PAACLwiVvAAAAAAwiUAIAAABgEIESAAAAAIMIlAAAAAAYRKAEAAAAwCC+5Q0AAIDBfBsobG3OUAIAAABgEIESAAAAAIMIlAAAAAAYRKAEAAAAwCACJQAAAAAGESgBAAAAMIhACQAAAIBBBEoAAAAADCJQAgAAAGAQgRIAAAAAgwiUAAAAABhEoAQAAADAIAIlAAAAAAYRKAEAAAAwiEAJAAAAgEEESgAAAAAMIlACAAAAYBCBEgAAAACDzG2gVFUPrKovV9WZVfUns64HAAAAgJG9Zl3ASqpqzySvSHK/JOcm+XRVvau7vzDbymAx3OcFb1vX7X3o2b++rtsDYLb0CQBgLXMZKCW5a5Izu/usJKmqf0ry0CQCJQAAAJghf3ggmd9A6cAk5yyZPzfJL86olp/yjQsuWtft3eDG11nX7QEwW/oEAJPoE8AiqO6edQ0/paoemeQB3f374/nHJ7lrd/+nJY85IskRSbJt27a7nH322UmSSy/6/rrWcq3rXHNdt7cVvfsjX1rX7T340Nuu6/bmzX981Qnrur3XPekB67q9eXLAE1++7ts8/5ij1n2b07je4U9e1+19+32vXDpb67rxTUKf2Dz0iemtd49I9IkhZtUjkg3tE1uyRyT6xGaiT0xPnxhGn1jdNH1iXs9QOjfJQUvmb5Lk/KUP6O6jkxydJNu3b/9xKuaADUCiTwAwmT4BsGvm9VvePp3kVlV186raO8ljkrxrxjUBAAAAkDk9Q6m7f1hVRyY5IcmeSV7b3WfMuCwAAAAAMqeBUpJ093uTvHfWdQAAADA7i3zPI2Zrlvc8WgRzGygBAAAjy75sAQCuYhZ9QqAEMAf8jwIAALCZCJQAAFhILmUAgI0zr9/yBgAAAMCccoYSsGn5yzMAAMBsCJQAAABgAbzuSQ+YdQlsIS55AwAAAGAQZygBAGwS/vIMAMwLZygBAAAAMIhACQAAAIBBXPIGc8blDAAAAMw7ZygBAAAAMIhACQAAAIBBBEoAAAAADCJQAgAAAGAQgRIAAAAAgwiUAAAAABhEoAQAAADAIAIlAAAAAAYRKAEAAAAwiEAJAAAAgEGmCpSq6kXTLAMAAABg8U17htL9Vlh2+HoWAgAAAMDmsNeklVX15CR/mOQWVXXaklX7Jfn4RhYGAAAAwHyaGCgleWOS9yV5QZI/WbL84u7+9oZVBQAAAMDcmhgodfd3k3w3yWOras8kNxz/zr5VtW93f3031AgAAADAHFnrDKUkSVUdmeTPkvxbkh+NF3eSO21MWQAAAADMq6kCpSRPTXKb7v7WRhYDAAAAwPyb9lvezsno0jcAAAAAtrhpz1A6K8mJVXVckst3LOzul25IVQAAAADMrWkDpa+Pp73HEwAAAABb1FSBUnf/eZJU1bW6+9KNLQkAAACAeTbVPZSq6u5V9YUkXxzP/3xV/d2GVgYAAADAXJr2ptwvS/KAJN9Kku4+NcmhG1UUAAAAAPNr2kAp3X3OskVXrnMtAAAAAGwC096U+5yq+qUkXVV7Jzkq48vfAAAAANhapj1D6UlJnpLkwCTnJjl4PA8AAADAFjPtt7xdmORxG1wLAAAAAJvAxECpqp7Z3S+uqv+RpJev7+6jNqwyAAAAAObSWmco7bhP0mc2uhAAAAAANoeJgVJ3v3v87z/unnIAAAAAmHdT3ZS7qj5QVddZMn/dqjph48oCAAAAYF5N+y1v1+/ui3bMdPd3ktxgY0oCAAAAYJ5NGyhdWVXbdsxU1U2zwk26AQAAAFh8a92Ue4fnJPlYVZ00nj80yREbUxIAAAAA82yqQKm7j6+qQ5LcLUkl+ePuvnBDKwMAAABgLk285K2qbjv+95Ak25Kcn+S8JNvGy9ZdVf1ZVZ1XVaeMpwdtxH4AAAAA2DlrnaH0tIwubfubFdZ1kvuse0Ujf9vdL9mgbQMAAACwC9YKlD4w/vcJ3X3WRhcDAAAAwPxb61venj3+960bXcgyR1bVaVX12qq67m7eNwAAAAATrHWG0rer6sNJblFV71q+srsfsjM7raoPJrnRCquek+SVSf4yo0vq/jKjy+1+b4VtHJHxN81t27ZtZ8oAYIHpEwBMok8A7Jq1AqUHJTkkyeuz8n2Udkp333eax1XVMUnes8o2jk5ydJJs376916s2ABaDPgHAJPoEwK5ZK1B6TXc/vqqO6e6TdkdBVXXj7r5gPPvwJKfvjv0CAAAAMJ21AqW7VNVNkzxufLZQLV3Z3d/egJpeXFUHZ3TJ29eS/MEG7AMAAACAnbRWoPSqJMcnuUWSk3PVQKnHy9dVdz9+vbcJAAAAwPqZ+C1v3f3y7r5dktd29y26++ZLpnUPkwAAAACYf2udoZQk6e4nV9U9k9yqu/+hqvZPsl93f3VjywMAFt2DD73trEsAAGCgiWco7VBVz0vyrCTPHi/aO8n/3qiiAAAAAJhfUwVKGX3b2kOSXJok3X1+kv02qigAAAAA5te0gdIV3d0Z3Yg7VXWtjSsJAAAAgHk2baD0lqp6dZLrVNUTk3wwyTEbVxYAAAAA82ram3K/pKrul+R7SW6T5E+7+wMbWhkAAAAAc2mqQGnstCRXH/986gbUAgAAAMAmMO23vD0qyaeSPDLJo5J8sqp+YyMLAwAAAGA+TXuG0nOS/EJ3fyNJqur6Gd1H6a0bVRgAAAAA82nam3LvsSNMGvvWgN8FAAAAYIFMe4bS8VV1QpI3jecfneS9G1MSAAAAAPNsYqBUVbdMcsPufkZVPSLJPZNUkk8kecNuqA8AAACAObPWZWsvS3JxknT327v7ad39xxmdnfSyjS4OAAAAgPmzVqB0s+4+bfnC7v5MkpttSEUAAAAAzLW1AqV9Jqy7xnoWAgAAAMDmsFag9OmqeuLyhVX1hCQnb0xJAAAAAMyztb7l7alJjq2qx+UnAdL2JHsnefhGFgYAAADAfJoYKHX3vyX5paq6d5I7jhcf190f2vDKAAAAAJhLa52hlCTp7g8n+fAG1wIAAADAJrDWPZQAAAAA4CoESgAAAAAMIlACAAAAYBCBEgAAAACDCJQAAAAAGESgBAAAAMAgAiUAAAAABhEoAQAAADCIQAkAAACAQQRKAAAAAAwiUAIAAABgEIESAAAAAIMIlAAAAAAYRKAEAAAAwCACJQAAAAAGESgBAAAAMIhACQAAAIBBBEoAAAAADCJQAgAAAGAQgRIAAAAAgwiUAAAAABhEoAQAAADAIAIlAAAAAAYRKAEAAAAwiEAJAAAAgEFmEihV1SOr6oyq+lFVbV+27tlVdWZVfbmqHjCL+gAAAABY3V4z2u/pSR6R5NVLF1bV7ZM8JskdkhyQ5INVdevuvnL3lwgAAADASmZyhlJ3f7G7v7zCqocm+afuvry7v5rkzCR33b3VAQAAADDJvN1D6cAk5yyZP3e8DAAAAIA5sWGXvFXVB5PcaIVVz+nud672ayss61W2f0SSI5Jk27ZtO1UjAItLnwBgEn0CYNdsWKDU3ffdiV87N8lBS+ZvkuT8VbZ/dJKjk2T79u0rhk4AbF36BACT6BMAu2beLnl7V5LHVNXVq+rmSW6V5FMzrgkAAACAJWYSKFXVw6vq3CR3T3JcVZ2QJN19RpK3JPlCkuOTPMU3vAEAAADMlw275G2S7j42ybGrrHt+kufv3ooAAAAAmNa8XfIGAAAAwJwTKAEAAAAwiEAJAAAAgEEESgAAAAAMIlACAAAAYBCBEgAAAACDCJQAAAAAGESgBAAAAMAgAiUAAAAABhEoAQAAADCIQAkAAACAQQRKAAAAAAwiUAIAAABgEIESAAAAAIMIlAAAAAAYRKAEAAAAwCACJQAAAAAGESgBAAAAMIhACQAAAIBBBEoAAAAADCJQAgAAAGAQgRIAAAAAgwiUAAAAABhEoAQAAADAIAIlAAAAAAYRKAEAAAAwiEAJAAAAgEEESgAAAAAMIlACAAAAYBCBEgAAAACDCJQAAAAAGESgBAAAAMAgAiUAAAAABhEoAQAAADCIQAkAAACAQQRKAAAAAAwiUAIAAABgEIESAAAAAIMIlAAAAAAYRKAEAAAAwCACJQAAAAAGESgBAAAAMIhACQAAAIBBBEoAAAAADCJQAgAAAGAQgRIAAAAAg8wkUKqqR1bVGVX1o6ravmT5zarqsqo6ZTy9ahb1AQAAALC6vWa039OTPCLJq1dY95XuPng31wMAAADAlGYSKHX3F5OkqmaxewAAAAB2wazOUJrk5lX1uSTfS/Lc7v7orAti1zz40NvOugQAAABgHW1YoFRVH0xyoxVWPae737nKr12QZFt3f6uq7pLkHVV1h+7+3grbPyLJEUmybdu29SobgAWhTwAwiT4BsGs2LFDq7vvuxO9cnuTy8c8nV9VXktw6yWdWeOzRSY5Oku3bt/euVQvAotEnAJhEnwDYNTP5lrfVVNX1q2rP8c+3SHKrJGfNtioAAAAAlppJoFRVD6+qc5PcPclxVXXCeNWhSU6rqlOTvDXJk7r727OoEQAAAICVzepb3o5NcuwKy9+W5G27vyIAAAAApjVXl7wBAAAAMP8ESgAAAAAMIlACAAAAYBCBEgAAAACDCJQAAAAAGKS6e9Y17LKq+maSswf+2v5JLtyAcna3RRlHsjhjWZRxJIszlkUZR7JzY7mwux+4EcVsFvrEQowjWZyxLMo4ksUZy6KMIxk+li3fIxJ9IosxjmRxxrIo40gWZyyLMo5kHfvEQgRKO6OqPtPd22ddx65alHEkizOWRRlHsjhjWZRxJIs1lnm3KM/1oowjWZyxLMo4ksUZy6KMI1msscy7RXmuF2UcyeKMZVHGkSzOWBZlHMn6jsUlbwAAAAAMIlACAAAAYJCtHCgdPesC1smijCNZnLEsyjiSxRnLoowjWayxzLtFea4XZRzJ4oxlUcaRLM5YFmUcyWKNZd4tynO9KONIFmcsizKOZHHGsijjSNZxLFv2HkoAAAAA7JytfIYSAAAAADth4QOlqjqoqr5aVdcbz193PH/Tqjq+qi6qqvfMus61rDGOk6vqlKo6o6qeNOta17LGWK4cj+WUqnrXrGudZMI4fmfJGE6pqh9U1cNmXe8ka7wmL6qq08fTo2dd63I78xmvqiOr6syq6qrafzaVX9VOjuM1VXVqVZ1WVW+tqn1nU/3mpk/Ml0XpEcni9InN3CMSfUKf2DWL0iMSfWIe6RPzQZ/Y+T6xJS55q6pnJrlldx9RVa9O8rXufkFV/UqSayb5g+7+tdlWubaVxpHkbzJ6HS8fv/inJ/ml7j5/hqWuacJrckl3b5r/2FltHEvWXy/JmUlu0t3fn1Wd01jl/XVakqcmOTzJ1ZOclOQ+3f29mRW6gqGf8aq6c5LvJDkxyfbuvnAWdS+3E+O49o7XoqpemuQb3f3CmRS/yekT82VRekSyOH1iM/eIRJ8Y/6xP7KRF6RGJPjGP9In5oE/sZJ/o7oWfklwtP3kzn5Fk7yXr7pXkPbOucVfHMV7/s0m+nuSAWde6s2NJcsmsa1vn1+SIJG+YdZ07O5Ykz0jy3CWPeU2SR8261iGvw6TPeEaNbv9Z178O46gkr0zyrFmPYbNO+sR8TYvSI6Z8TTZFn9jMPWKt10GfMG3U8z6Pkz4xf5M+MR+TPrFzfWKvbAHd/e9V9Ywkxye5f3dfMeuadsZq46iqg5Icl+SWSZ7Rc/zXhB0mvCb7VNVnkvwwyQu7+x0zK3IKU7y3HpPkpbu/suFWGktVnZrkeeO0+ppJ7p3kC7OscyWL/hmfpKr+IcmDMnpd/vMGl7iwFv09tNn6xKL0iGRx+sRm7hHJ4n/GJ9Endt2ivH8SfWIe6RPzYVE+57u7Tyz8PZSWODzJBUnuOOtCdtFPjaO7z+nuO2XUAH6nqm44q+IGWuk12dbd25P8ZpKXVdXPzaSyYVZ8b1XVjZP8hyQnzKKonXSVsXT3+5O8N8n/TfKmJJ/IqEHPo4X9jE/S3b+b5IAkX0wyl9elbyIL+x7apH1iUXpEsjh9YjP3iGSBP+OT6BPrZlHeP4k+MY/0ifmwKJ/z3dYntkSgVFUHJ7lfkrsl+ePxB3PTWWsc478knJHkl2dQ3iCrjWXHX0O6+6yMrke986xqnMYar8mjkhzb3f8+k+IGmvCaPL+7D+7u+2V0KuT/m2GZK9oqn/HVdPeVSd6c5Nc3sLyFtlXeQ5ulTyxKj0gWp09s5h6RbJ3P+Gr0iV2zKO+fRJ+YR/rEfFiUz/nu7hMLHyhV1Y5rAZ/a3V9P8t+SvGS2VQ232jiq6iZVdY3xY66b5B5Jvjy7Stc2YSzXraqrjx+zf0ZjmctTIpOp3luPzSiJn3sTXpM9q+pnx4+5U5I7JXn/7Cr9aYv+GZ/0+Kq65ZLffXCSL+2OWhfNor+HNlufWJQekSxOn9jMPSJZ/M/4pMfrE7tuUd4/iT4xj/SJ+bAon/OZ9IkhN1zajFNGNzF785L5PZOcnOSwJB9N8s0klyU5N8kDZl3vTozjeRnddOvU8b9HzLrWXXxNPj8ey+eTPGHWte7COG6W5Lwke8y6znUYyxfG0z8nOXjWtQ6sfcXPeJKjxvM/THJ+kr/fbOPI6A8CHx9/Vk5P8oYk1571ODbjpE/M17QoPWKKsWyaPrGZe8QU9esTpnV//8zrpE/M36RPzMekT+x8n6jxjgAAAABgKgt/yRsAAAAA60ugBAAAAMAgAiUAAAAABhEoAQAAADCIQAkAAACAQfaadQGwu1XVlRl9NeJeSb6a5PHdfdFsqwJgXugTAEyiT8CIM5TYii7r7oO7+45Jvp3kKeux0arabQFtVe25u/YFsAXpEwBMok9ABErwiSQH7pipqmdU1aer6rSq+vMly/9rVX2pqj5QVW+qqqePl59YVX9dVScl+aOqun5VvW28jU9X1T3Gjzusqk4ZT5+rqv2q6sZV9ZHxstOr6pfHj31sVX1+vOxFS2q4pKr+oqo+meTuu+n5Adjq9AkAJtEn2LJc8saWNU7lfyXJa8bz909yqyR3TVJJ3lVVhyb5fpJfT3LnjD4zn01y8pJNXae7Dxtv441J/ra7P1ZV25KckOR2SZ6e5Cnd/fGq2jfJD5IckeSE7n7+uJZrVtUBSV6U5C5JvpPk/VX1sO5+R5JrJTm9u/90454VAHbQJwCYRJ9gqxMosRVdo6pOSXKzjA7kHxgvv/94+tx4ft+MGsJ+Sd7Z3ZclSVW9e9n23rzk5yzPwmgAAAH3SURBVPsmuX1V7Zi/dlXtl+TjSV5aVW9I8vbuPreqPp3ktVV1tSTv6O5Tquo+SU7s7m+O9/WGJIcmeUeSK5O8bT2eAAAm0icAmESfgLjkja3psu4+OMlNk+ydn1zzXEleML4e+uDuvmV3v2a8fJJLl/y8R5K7L9nGgd19cXe/MMnvJ7lGkn+uqtt290cyOrifl+T1VfXba+zrB9195eDRAjCUPgHAJPoERKDEFtbd301yVJKnj1P9E5L83vgU0lTVgVV1gyQfS/LgqtpnvO5XJ2z2/UmO3DFTVQeP//257v58d78oyWeS3LaqbprkG919TEanyR6S5JNJDquq/cenrT42yUnrO3IApqFPADCJPsFW55I3trTu/lxVnZrkMd39+qq6XZJPjE8xvSTJb3X3p6vqXUlOTXJ2Rgfw766yyaOSvKKqTsvo8/WRJE9K8tSqundGp5l+Icn7kjwmyTOq6t/H+/rt7r6gqp6d5MMZ/XXhvd39zg0ZPABr0icAmESfYCur7p51DTD3qmrf7r6kqq6Z0UH9iO7+7KzrAmA+6BMATKJPsIicoQTTObqqbp9knyT/6OAPwDL6BACT6BMsHGcoAQAAADCIm3IDAAAAMIhACQAAAIBBBEoAAAAADCJQAgAAAGAQgRIAAAAAgwiUAAAAABjk/wNVFd+s0+VaaQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x324 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Xs = ['X'+str(i + 1) for i in range(X.shape[1])]\n",
    "taus = [100, 10, 1]\n",
    "\n",
    "fig, ax = plt.subplots(ncols = len(taus), figsize = (20, 4.5), sharey = True)\n",
    "for i, tau in enumerate(taus):\n",
    "    model = BayesianRegression()\n",
    "    model.fit(X, y, sigma_squared, tau) \n",
    "    betas = model.beta_hats[1:]\n",
    "    sns.barplot(Xs, betas, ax = ax[i], palette = 'PuBu')\n",
    "    ax[i].set(xlabel = 'Regressor', title = fr'Regression Coefficients with $\\tau = $ {tau}')\n",
    "    ax[i].set(xticks = np.arange(0, len(Xs), 2), xticklabels = Xs[::2])\n",
    "\n",
    "ax[0].set(ylabel = 'Coefficient')\n",
    "sns.set_context(\"talk\")\n",
    "sns.despine();"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
